<!DOCTYPE html>
<html class="writer-html5" lang="en" >
<head>
    <meta charset="utf-8" />
    <meta http-equiv="X-UA-Compatible" content="IE=edge" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="author" content="Patrick T. Komiske III" /><link rel="canonical" href="https://thaler-lab.github.io/Wasserstein/docs/emd/" />
      <link rel="shortcut icon" href="../../img/favicon.ico" />
    <title>EMD - Wasserstein</title>
    <link rel="stylesheet" href="../../css/theme.css" />
    <link rel="stylesheet" href="../../css/theme_extra.css" />
        <link href="../../css/eftheme.css" rel="stylesheet" />
        <link href="../../css/pygmentize_friendly.css" rel="stylesheet" />

      <script>
        // Current page data
        var mkdocs_page_name = "EMD";
        var mkdocs_page_input_path = "docs/emd.md";
        var mkdocs_page_url = "/Wasserstein/docs/emd/";
      </script>

    <script src="../../js/jquery-3.6.0.min.js" defer></script>
    <!--[if lt IE 9]>
      <script src="../../js/html5shiv.min.js"></script>
    <![endif]-->
      <script>
        (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
        (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
        m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
        })(window,document,'script','https://www.google-analytics.com/analytics.js','ga');

        ga('create', 'UA-122962541-3', 'Wasserstein');
        ga('send', 'pageview');
      </script>

</head>

<body class="wy-body-for-nav" role="document">

  <div class="wy-grid-for-nav">
    <nav data-toggle="wy-nav-shift" class="wy-nav-side stickynav">
    <div class="wy-side-scroll">
      <div class="wy-side-nav-search">
  <a href="../..">
    <div class="eflogo">
      <img src="../../img/eflogowhite.png"  class="eflogo-img"> Wasserstein
    </div>
  </a><div role="search">
  <form id ="rtd-search-form" class="wy-form" action="../../search.html" method="get">
      <input type="text" name="q" placeholder="Search docs" title="Type search term here" />
  </form>
</div>
      </div>

      <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="Navigation menu">
              <ul>
                <li class="toctree-l1"><a class="reference internal" href="../..">Home</a>
                </li>
              </ul>
              <p class="caption"><span class="caption-text">Getting Started</span></p>
              <ul>
                  <li class="toctree-l1"><a class="reference internal" href="../../installation/">Installation</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../../demos/">Python Demos</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../../examples/">C++ Examples</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../../faqs/">FAQs</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../../releases/">Release Notes</a>
                  </li>
              </ul>
              <p class="caption"><span class="caption-text">Documentation</span></p>
              <ul class="current">
                  <li class="toctree-l1 current"><a class="reference internal current" href="./">EMD</a>
    <ul class="current">
    <li class="toctree-l2"><a class="reference internal" href="#python">Python</a>
        <ul>
    <li class="toctree-l3"><a class="reference internal" href="#run-computation">Run Computation</a>
        <ul>
    <li class="toctree-l4"><a class="reference internal" href="#__call__">__call__</a>
    </li>
        </ul>
    </li>
    <li class="toctree-l3"><a class="reference internal" href="#access-results">Access Results</a>
        <ul>
    <li class="toctree-l4"><a class="reference internal" href="#emd_1">emd</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#dists">dists</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#flows">flows</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#flow">flow</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#status">status</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#n0">n0</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#n1">n1</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#extra">extra</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#weightdiff">weightdiff</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#scale">scale</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#duration">duration</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#n_iter">n_iter</a>
    </li>
        </ul>
    </li>
    <li class="toctree-l3"><a class="reference internal" href="#getset-options">Get/Set Options</a>
    </li>
    <li class="toctree-l3"><a class="reference internal" href="#other-methods">Other Methods</a>
        <ul>
    <li class="toctree-l4"><a class="reference internal" href="#preprocess_centerweightedcentroid">preprocess_CenterWeightedCentroid</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#description">description</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#clear">clear</a>
    </li>
        </ul>
    </li>
        </ul>
    </li>
    <li class="toctree-l2"><a class="reference internal" href="#c">C++</a>
        <ul>
    <li class="toctree-l3"><a class="reference internal" href="#run-computation_1">Run Computation</a>
        <ul>
    <li class="toctree-l4"><a class="reference internal" href="#operator">operator()</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#compute">compute</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#ground_dists">ground_dists</a>
    </li>
        </ul>
    </li>
    <li class="toctree-l3"><a class="reference internal" href="#access-results_1">Access Results</a>
    </li>
    <li class="toctree-l3"><a class="reference internal" href="#getset-options_1">Get/Set Options</a>
    </li>
    <li class="toctree-l3"><a class="reference internal" href="#other-methods_1">Other Methods</a>
        <ul>
    <li class="toctree-l4"><a class="reference internal" href="#preprocess">preprocess</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#description_1">description</a>
    </li>
    <li class="toctree-l4"><a class="reference internal" href="#clear_1">clear</a>
    </li>
        </ul>
    </li>
        </ul>
    </li>
    </ul>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../emds/">PairwiseEMD</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../externalemdhandler/">External EMD Handlers</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../event/">Events</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../pairwisedistance/">Pairwise Distance</a>
                  </li>
                  <li class="toctree-l1"><a class="reference internal" href="../utils/">Utils</a>
                  </li>
              </ul>
      </div>
    </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">
      <nav class="wy-nav-top" role="navigation" aria-label="Mobile navigation menu">
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="../..">Wasserstein</a>

      </nav>
      <div class="wy-nav-content">
        <div class="rst-content"><div role="navigation" aria-label="breadcrumbs navigation">
  <ul class="wy-breadcrumbs">
    <li><a href="../.." class="icon icon-home" alt="Docs"></a> &raquo;</li>
          <li>Documentation &raquo;</li><li>EMD</li>
    <li class="wy-breadcrumbs-aside">
    </li>
  </ul>
  <hr/>
</div>

          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
            <div class="section" itemprop="articleBody">

                <h1 id="emd">EMD</h1>
<p>The <code>EMD</code> class handles computation of individual Wasserstein distances between pairs of distributions. Adopting language from particle physics, we will call the distributions "events," the discrete entities in the ground space "particles," and the particle weights (probability mass) "energy".</p>
<p>The optimal transport problem can be specified as:</p>
<div>
<div class="MathJax_Preview">\text{EMD}_{\beta,R}(\mathcal E,\mathcal E^\prime)=\min_{\{f_{ij}\ge0\}}\sum_{i=1}^M\sum_{j=1}^{M'}f_{ij}\left(\frac{
\theta_{ij}}{R}\right)^\beta + \left|\sum_{i=1}^ME_i-\sum_{j=1}^{M'}E^\prime_j
\right|,</div>
<script type="math/tex; mode=display">\text{EMD}_{\beta,R}(\mathcal E,\mathcal E^\prime)=\min_{\{f_{ij}\ge0\}}\sum_{i=1}^M\sum_{j=1}^{M'}f_{ij}\left(\frac{
\theta_{ij}}{R}\right)^\beta + \left|\sum_{i=1}^ME_i-\sum_{j=1}^{M'}E^\prime_j
\right|,</script>
</div>
<div>
<div class="MathJax_Preview">\sum_{j=1}^{M'}f_{ij}\le E_i, \quad \sum_{i=1}^Mf_{ij}\le E^\prime_j,
\quad\sum_{i=1}^M\sum_{j=1}^{M'}f_{ij}=\min\left(\sum_{i=1}^ME_i,\sum_{j=1}^{M'}E^\prime_j\right),</div>
<script type="math/tex; mode=display">\sum_{j=1}^{M'}f_{ij}\le E_i, \quad \sum_{i=1}^Mf_{ij}\le E^\prime_j,
\quad\sum_{i=1}^M\sum_{j=1}^{M'}f_{ij}=\min\left(\sum_{i=1}^ME_i,\sum_{j=1}^{M'}E^\prime_j\right),</script>
</div>
<p>where <span><span class="MathJax_Preview">\mathcal E</span><script type="math/tex">\mathcal E</script></span> and <span><span class="MathJax_Preview">\mathcal E'</span><script type="math/tex">\mathcal E'</script></span> are events with <span><span class="MathJax_Preview">M</span><script type="math/tex">M</script></span> and <span><span class="MathJax_Preview">M'</span><script type="math/tex">M'</script></span> particles, energies <span><span class="MathJax_Preview">\{E_i\}</span><script type="math/tex">\{E_i\}</script></span> and <span><span class="MathJax_Preview">\{E'_j\}</span><script type="math/tex">\{E'_j\}</script></span>, respectively, and pairwise distance matrix <span><span class="MathJax_Preview">\theta_{ij}</span><script type="math/tex">\theta_{ij}</script></span>. <span><span class="MathJax_Preview">R&gt;0</span><script type="math/tex">R>0</script></span> is a parameter controlling the tradeoff between transporting and creating/destroying energy (<span><span class="MathJax_Preview">R</span><script type="math/tex">R</script></span> should be large than half the maximum distance in the ground space in order to ensure the EMD defines a proper metric). <span><span class="MathJax_Preview">\beta&gt;0</span><script type="math/tex">\beta>0</script></span> is an angular exponent, and the EMD only satisfies the triangle inequality for <span><span class="MathJax_Preview">\beta&gt;1</span><script type="math/tex">\beta>1</script></span> if it is raised to the <span><span class="MathJax_Preview">1/\beta</span><script type="math/tex">1/\beta</script></span> power.</p>
<hr />
<h2 id="python">Python</h2>
<p>The Python <code>EMD</code> function returns an object (either <code>wasserstein.EMDFloat64</code> or <code>wasserstein.EMDFloat32</code>, each of which are instantiations of the C++ template <code>EMD</code>) that can be used to compute EMD distances. It is designed to work with numpy arrays efficiently. The <code>EMDYPhi</code> function behaves similarly but implements <span><span class="MathJax_Preview">2\pi</span><script type="math/tex">2\pi</script></span> periodicity in the second coordinate dimension, and so is suited for using the rapidity-azimuth plane as a ground space.</p>
<div class="codehilite"><pre><span></span><code><span class="n">wasserstein</span><span class="o">.</span><span class="n">EMD</span><span class="p">(</span><span class="n">R</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">beta</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">norm</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span> <span class="n">do_timing</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span>
                <span class="n">n_iter_max</span><span class="o">=</span><span class="mi">100000</span><span class="p">,</span> <span class="n">epsilon_large_factor</span><span class="o">=</span><span class="mf">1000.0</span><span class="p">,</span> <span class="n">epsilon_small_factor</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span>
                <span class="n">dtype</span><span class="o">=</span><span class="s1">&#39;float64&#39;</span><span class="p">)</span>

<span class="n">wasserstein</span><span class="o">.</span><span class="n">EMDYPhi</span><span class="p">(</span><span class="n">R</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">beta</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">norm</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span> <span class="n">do_timing</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span>
                    <span class="n">n_iter_max</span><span class="o">=</span><span class="mi">100000</span><span class="p">,</span> <span class="n">epsilon_large_factor</span><span class="o">=</span><span class="mf">1000.0</span><span class="p">,</span> <span class="n">epsilon_small_factor</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span>
                    <span class="n">dtype</span><span class="o">=</span><span class="s1">&#39;float64&#39;</span><span class="p">)</span>
</code></pre></div>

<p><strong>Arguments</strong></p>
<ul>
<li><strong>R</strong> : <em>float</em><ul>
<li>The R parameter in the EMD definition that controls the relative
importance of the two terms. Must be greater than or equal to half
of the maximum ground distance in the space in order for the EMD
to be a valid metric satisfying the triangle inequality.</li>
</ul>
</li>
<li><strong>beta</strong> : <em>float</em><ul>
<li>The angular weighting exponent. The internal pairwsie distance
matrix is raised to this power prior to solving the optimal
transport problem.</li>
</ul>
</li>
<li><strong>norm</strong> : <em>bool</em><ul>
<li>Whether or not to normalize the particle weights to sum to one
prior to computing the EMD.</li>
</ul>
</li>
<li><strong>do_timing</strong> : <em>bool</em><ul>
<li>Whether or not to keep track of the duration of the underlying computation.</li>
</ul>
</li>
<li><strong>n_iter_max</strong> : <em>int</em><ul>
<li>Maximum number of iterations for solving the optimal transport
problem.</li>
</ul>
</li>
<li><strong>epsilon_large_factor</strong> : <em>float</em><ul>
<li>Controls some tolerances in the optimal transport solver. This
value is multiplied by the floating points epsilon (around 1e-16 for
64-bit floats) to determine the actual tolerance.</li>
</ul>
</li>
<li><strong>epsilon_small_factor</strong> : <em>float</em><ul>
<li>Analogous to <code>epsilon_large_factor</code> but used where the numerical
tolerance can be stricter.</li>
</ul>
</li>
<li><strong>dtype</strong> : one of <code>{'float64', 'float32', numpy.float32, numpy.float64}</code><ul>
<li>Controls the floating point precision used in the EMD computation.</li>
</ul>
</li>
</ul>
<hr />
<h3 id="run-computation">Run Computation</h3>
<h4 id="__call__"><strong><em>_</em>_call</strong>__</h4>
<div class="codehilite"><pre><span></span><code><span class="fm">__call__</span><span class="p">(</span><span class="n">weights0</span><span class="p">,</span> <span class="n">coords0</span><span class="p">,</span> <span class="n">weights1</span><span class="p">,</span> <span class="n">coords1</span><span class="p">)</span>
<span class="fm">__call__</span><span class="p">(</span><span class="n">weights0</span><span class="p">,</span> <span class="n">weights1</span><span class="p">,</span> <span class="n">dists</span><span class="p">)</span>
</code></pre></div>

<p>Compute the Wasserstein distance between two events using the EMD object.</p>
<p><strong>Arguments</strong></p>
<ul>
<li><strong>weights0</strong> : <em>1d numpy.ndarray</em><ul>
<li>The weight values of the first event.</li>
</ul>
</li>
<li><strong>weights1</strong> : <em>1d numpy.ndarray</em><ul>
<li>The weight values of the second event.</li>
</ul>
</li>
</ul>
<p>In the first version, the Euclidean ground distance is used where the particle coordinates are specified as:</p>
<ul>
<li><strong>coords0</strong> : <em>2d numpy.ndarray</em><ul>
<li>The Cartesian coordinates of the particles in the first event. If the ground space dimension is <code>d</code>, then it has shape <code>(M0, d)</code> where <code>M0</code> is the length of <code>weights0</code>.</li>
</ul>
</li>
<li><strong>coords1</strong> : <em>2d numpy.ndarray</em><ul>
<li>The Cartesian coordinates of the particles in the second event. If the ground space dimension is <code>d</code>, then it has shape <code>(M1, d)</code> where <code>M1</code> is the length of <code>weights1</code>.</li>
</ul>
</li>
</ul>
<p>In the second version, the ground space distances are provided directly:</p>
<ul>
<li><strong>dists</strong> : <em>2d numpy.ndarray</em><ul>
<li>The ground distances between particles in the first and second events. To use the above notation, it has shape <code>(M0, M1)</code>.</li>
</ul>
</li>
</ul>
<p><strong>Returns</strong></p>
<ul>
<li><em>float</em><ul>
<li>The Wasserstein distance between event0 and event1. See also <code>emd()</code>.</li>
</ul>
</li>
</ul>
<hr />
<h3 id="access-results">Access Results</h3>
<p>These methods access quantities determined during the most recent computation.</p>
<h4 id="emd_1">emd</h4>
<div class="codehilite"><pre><span></span><code><span class="n">emd</span><span class="p">()</span>
</code></pre></div>

<p>The Wasserstein distance.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>float</em><ul>
<li>The cost of transporting <code>event0</code> to <code>event1</code>.</li>
</ul>
</li>
</ul>
<h4 id="dists">dists</h4>
<div class="codehilite"><pre><span></span><code><span class="n">dists</span><span class="p">()</span>
</code></pre></div>

<p>The ground distances between the particles.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>numpy.ndarray</em> (first version only)<ul>
<li>The ground distance matrix as an array with shape <code>(n0, n1)</code>.</li>
</ul>
</li>
</ul>
<h4 id="flows">flows</h4>
<div class="codehilite"><pre><span></span><code><span class="n">flows</span><span class="p">()</span>
</code></pre></div>

<p>The optimal flow matrix.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>numpy.ndarray</em> (first version only)<ul>
<li>The flow matrix as an array with shape <code>(n0, n1)</code>.</li>
</ul>
</li>
</ul>
<h4 id="flow">flow</h4>
<div class="codehilite"><pre><span></span><code><span class="n">flow</span><span class="p">(</span><span class="n">i</span><span class="p">,</span> <span class="n">j</span><span class="p">)</span>
<span class="n">flow</span><span class="p">(</span><span class="n">ind</span><span class="p">)</span>
</code></pre></div>

<ul>
<li>First version: Access the flow between particle <code>i</code> in <code>event0</code> and <code>j</code> in <code>event1</code>. Note that negative indices are accepted and count from the end, as usual in Python.</li>
<li>Second version: Access the flow using a 1D index such that <code>ind = i*n1() + j</code>.</li>
</ul>
<p><strong>Arguments</strong></p>
<ul>
<li><strong>i</strong> : <em>int</em><ul>
<li>Index of particle in <code>event0</code>. An <code>IndexError</code> is raised if out of range.</li>
</ul>
</li>
<li><strong>j</strong> : <em>int</em><ul>
<li>Index of particle in <code>event1</code>. An <code>IndexError</code> is raised if out of range.</li>
</ul>
</li>
<li><strong>ind</strong> : <em>int</em><ul>
<li>Raw flow index, equal to <code>i*n1() + j</code>.</li>
</ul>
</li>
</ul>
<p><strong>Returns</strong></p>
<ul>
<li><em>float</em><ul>
<li>The amount of flow between the specified particles.</li>
</ul>
</li>
</ul>
<h4 id="status">status</h4>
<div class="codehilite"><pre><span></span><code><span class="n">status</span><span class="p">()</span>
</code></pre></div>

<p>The status code of the solver. A non-zero value indicates a problem. See the <code>check_emd_status</code> function for interpretting the code.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>int</em><ul>
<li>Networks simplex status code. Zero indicates success, a non-zero value indicates an error.</li>
</ul>
</li>
</ul>
<h4 id="n0">n0</h4>
<div class="codehilite"><pre><span></span><code><span class="n">n0</span><span class="p">()</span>
</code></pre></div>

<p>Number of entities passed to the network simplex solver for <code>event0</code>.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>int</em><ul>
<li>The number of particles ultimately used in <code>event0</code>. This may be one greater than the number of given particles in <code>event0</code> if <code>norm=False</code> and <code>event0</code> has less total weight than <code>event1</code>.</li>
</ul>
</li>
</ul>
<h4 id="n1">n1</h4>
<div class="codehilite"><pre><span></span><code><span class="n">n1</span><span class="p">()</span>
</code></pre></div>

<p>Number of entities passed to the network simplex solver for <code>event1</code>.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>int</em><ul>
<li>The number of particles ultimately used in <code>event1</code>. This may be one greater than the number of given particles in <code>event1</code> if <code>norm=False</code> and <code>event1</code> has less total weight than <code>event0</code>.</li>
</ul>
</li>
</ul>
<h4 id="extra">extra</h4>
<div class="codehilite"><pre><span></span><code><span class="n">extra</span><span class="p">()</span>
</code></pre></div>

<p>Which event, if any, got an extra particle.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>int</em><ul>
<li>Which event, <code>0</code> or <code>1</code>, got an extra particle, or if neither did, <code>-1</code>.</li>
</ul>
</li>
</ul>
<h4 id="weightdiff">weightdiff</h4>
<div class="codehilite"><pre><span></span><code><span class="n">weightdiff</span><span class="p">()</span>
</code></pre></div>

<p>Total weight in <code>event1</code> minus total weight in <code>event0</code>.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>float</em><ul>
<li>The internally calculated difference in total weight between <code>event1</code> and <code>event0</code>.</li>
</ul>
</li>
</ul>
<h4 id="scale">scale</h4>
<div class="codehilite"><pre><span></span><code><span class="n">scale</span><span class="p">()</span>
</code></pre></div>

<p>The optimal transport problem is internally rescaled in order to provide more numerically stable computations. The <code>scale</code> is the value that the weights are divided by prior to running the network simplex algorithm.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>float</em><ul>
<li>The internally utilized scaling factor for the weights.</li>
</ul>
</li>
</ul>
<h4 id="duration">duration</h4>
<div class="codehilite"><pre><span></span><code><span class="n">duration</span><span class="p">()</span>
</code></pre></div>

<p>The duration of the core of the optimal transport calculation via the network simplex algorithm.</p>
<p><strong>Returns</strong></p>
<ul>
<li><em>float</em><ul>
<li>The time, in seconds, of the computation. If <code>do_timing</code> was not <code>True</code> for this computation then this value is not meaningful.</li>
</ul>
</li>
</ul>
<h4 id="n_iter">n_iter</h4>
<p>````python
n_iter()</p>
<div class="codehilite"><pre><span></span><code>The number of iterations the most recent optimal transport calculation took to finish.

**Returns**

- _int_
    - Number of iterations run by the Network Simplex solver.

#### node_potentials

```python
node_potentials()
</code></pre></div>

<p>The potentials assigned in the Network Simplex algorithm to the particles in each event.</p>
<p><strong>Returns</strong></p>
<ul>
<li>(<em>numpy.ndarray</em>, <em>numpy.ndarray</em>)<ul>
<li>A pair of NumPy arrays corresponding to the node potentials in the first and second event.</li>
</ul>
</li>
</ul>
<hr />
<h3 id="getset-options">Get/Set Options</h3>
<p>These methods can be used to get the current settings of the <code>EMD</code> object or to set new ones. The getter methods are:</p>
<ul>
<li><code>R()</code> : returns <em>float</em></li>
<li><code>beta()</code> : returns <em>float</em></li>
<li><code>norm()</code> : returns <em>bool</em></li>
<li><code>do_timing</code> : returns <em>bool</em></li>
</ul>
<p>The setter methods are:</p>
<ul>
<li><code>set_R(new_R)</code> : accepts <em>float</em></li>
<li><code>set_beta(new_beta)</code> : accepts <em>float</em></li>
<li><code>set_norm(new_norm)</code> : accepts <em>bool</em></li>
<li><code>set_do_timing(new_do_timing)</code> : accepts <em>bool</em></li>
<li><code>set_network_simplex_params(n_iter_max=100000, epsilon_large_factor=1000.0, epsilon_small_factor=1.0)</code><ul>
<li>This method resets all of the underlying network simplex solver's parameters at once.</li>
</ul>
</li>
</ul>
<hr />
<h3 id="other-methods">Other Methods</h3>
<h4 id="preprocess_centerweightedcentroid">preprocess_CenterWeightedCentroid</h4>
<div class="codehilite"><pre><span></span><code><span class="n">preprocess_CenterWeightedCentroid</span><span class="p">()</span>
</code></pre></div>

<p>This method adds a <code>CenterWeightedCentroid</code> preprocessor to the internal list. The particles of each event will be adjusted so that the origin of the ground space corresponds to their weighted centroid.</p>
<h4 id="description">description</h4>
<div class="codehilite"><pre><span></span><code><span class="n">description</span><span class="p">(</span><span class="n">write_preprocessors</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
</code></pre></div>

<p>Returns a string that describes the <code>EMD</code> object. Printing the <code>EMD</code> object uses this method to describe the object.</p>
<p><strong>Arguments</strong></p>
<ul>
<li><strong>write_preprocessors</strong> : <em>bool</em><ul>
<li>Whether or not to include preprocessors in the description. There are currently no preprocessors included in the Wasserstein package but this may change in the future.</li>
</ul>
</li>
</ul>
<p><strong>Returns</strong></p>
<ul>
<li><em>string</em><ul>
<li>The description of the <code>EMD</code> object.</li>
</ul>
</li>
</ul>
<h4 id="clear">clear</h4>
<div class="codehilite"><pre><span></span><code><span class="n">clear</span><span class="p">()</span>
</code></pre></div>

<p>Frees some memory in use by the <code>EMD</code> object. This should not normally need to be called by the user.</p>
<hr />
<h2 id="c">C++</h2>
<p>The Wasserstein C++ interface makes heavy use of templates in order to give the user full control over the EMD computation. Most classes are templated to accept a floating point type (either <code>float</code> or <code>double</code>), though some aliases are provided to make things easier. The <code>EMD</code> class provides the implementation of the EMD as defined above. It accepts four template parameters: a floating point <code>Value</code> type, an <code>Event</code> class, a <code>PairwiseDistance</code> class, and a <code>NetworkSimplex</code> class. The default namespace for the package is <code>wasserstein</code>.</p>
<div class="codehilite"><pre><span></span><code><span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span> <span class="nc">Value</span><span class="p">,</span>
         <span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">Event</span> <span class="o">=</span> <span class="n">DefaultEvent</span><span class="p">,</span>
         <span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">PairwiseDistance</span> <span class="o">=</span> <span class="n">DefaultPairwiseDistance</span><span class="p">,</span>
         <span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">NetworkSimplex</span> <span class="o">=</span> <span class="n">DefaultNetworkSimplex</span><span class="o">&gt;</span>
<span class="n">wasserstein</span><span class="o">::</span><span class="n">EMD</span><span class="p">(</span><span class="n">Value</span> <span class="n">R</span> <span class="o">=</span> <span class="mi">1</span><span class="p">,</span> <span class="n">Value</span> <span class="n">beta</span> <span class="o">=</span> <span class="mi">1</span><span class="p">,</span> <span class="kt">bool</span> <span class="n">norm</span> <span class="o">=</span> <span class="nb">false</span><span class="p">,</span>
                 <span class="kt">bool</span> <span class="n">do_timing</span> <span class="o">=</span> <span class="nb">false</span><span class="p">,</span> <span class="kt">bool</span> <span class="n">external_dists</span> <span class="o">=</span> <span class="nb">false</span><span class="p">,</span>
                 <span class="kt">unsigned</span> <span class="n">n_iter_max</span> <span class="o">=</span> <span class="mi">100000</span><span class="p">,</span>
                 <span class="n">Value</span> <span class="n">epsilon_large_factor</span> <span class="o">=</span> <span class="mi">1000</span><span class="p">,</span>
                 <span class="n">Value</span> <span class="n">epsilon_small_factor</span> <span class="o">=</span> <span class="mi">1</span><span class="p">);</span>
</code></pre></div>

<p>See the <a href="#python">Python class</a> of the same name for the meaning of these arguments. <code>external_dists</code> is a new flag indicating whether or not the internal <code>PairwiseDistance</code> object should be used to obtain ground distances or whether the ground distances will be externally provided (see the <code>ground_dists()</code> method below).</p>
<p>Ordinarily, users will make use of one of the template specializations of <code>EMD</code>, which only require one to specify an <code>Event</code> and a <code>PairwiseDistance</code> class. These are:</p>
<div class="codehilite"><pre><span></span><code><span class="c1">// double precision, default NetworkSimplex</span>
<span class="k">template</span><span class="o">&lt;</span><span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">Event</span> <span class="o">=</span> <span class="n">DefaultEvent</span><span class="p">,</span>
         <span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">PairwiseDistance</span> <span class="o">=</span> <span class="n">DefaultPairwiseDistance</span><span class="o">&gt;</span>
<span class="k">using</span> <span class="n">EMDFloat64</span> <span class="o">=</span> <span class="n">EMD</span><span class="o">&lt;</span><span class="kt">double</span><span class="p">,</span> <span class="n">Event</span><span class="p">,</span> <span class="n">PairwiseDistance</span><span class="o">&gt;</span><span class="p">;</span>

<span class="c1">// single precision, default NetworkSimplex</span>
<span class="k">template</span><span class="o">&lt;</span><span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">Event</span> <span class="o">=</span> <span class="n">DefaultEvent</span><span class="p">,</span>
         <span class="k">template</span><span class="o">&lt;</span><span class="k">typename</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">PairwiseDistance</span> <span class="o">=</span> <span class="n">DefaultPairwiseDistance</span><span class="o">&gt;</span>
<span class="k">using</span> <span class="n">EMDFloat32</span> <span class="o">=</span> <span class="n">EMD</span><span class="o">&lt;</span><span class="kt">float</span><span class="p">,</span> <span class="n">Event</span><span class="p">,</span> <span class="n">PairwiseDistance</span><span class="o">&gt;</span><span class="p">;</span>
</code></pre></div>

<p>See the <a href="../event">Event</a> and <a href="../pairwisedistance">PairwiseDistance</a> docs for more on what classes are provided to fill these roles.</p>
<hr />
<h3 id="run-computation_1">Run Computation</h3>
<p>There are two methods that can be used to compute an EMD distance between two events, <code>operator()</code> and <code>compute</code>. The difference is that the first is templated to accept "proto events" (which may be of type <code>Event</code>, or are something that <code>Event</code> can be constructed from) and also preprocesses the events. The second accepts two already constructed events, does not preprocess them, and does not check that their weights have been setup properly.</p>
<h4 id="operator">operator()</h4>
<div class="codehilite"><pre><span></span><code><span class="k">template</span><span class="o">&lt;</span><span class="k">class</span> <span class="nc">ProtoEvent0</span><span class="p">,</span> <span class="k">class</span> <span class="nc">ProtoEvent1</span><span class="o">&gt;</span>
<span class="n">Value</span> <span class="k">operator</span><span class="p">()(</span><span class="k">const</span> <span class="n">ProtoEvent0</span> <span class="o">&amp;</span> <span class="n">pev0</span><span class="p">,</span> <span class="k">const</span> <span class="n">ProtoEvent1</span> <span class="o">&amp;</span> <span class="n">pev1</span><span class="p">);</span>
</code></pre></div>

<p>This version preprocesses each event to ensure any preprocessors are called and the weights are normalized properly (if <code>norm=true</code>).</p>
<p><strong>Arguments</strong></p>
<ul>
<li><strong>pev0</strong><ul>
<li>The first event will be constructed by calling <code>Event(pev0)</code>.</li>
</ul>
</li>
<li><strong>pev1</strong><ul>
<li>The second event will be constructed by calling <code>Event(pev1)</code>.</li>
</ul>
</li>
</ul>
<p><strong>Returns</strong></p>
<ul>
<li>The Wasserstein distance between the two events.</li>
</ul>
<h4 id="compute">compute</h4>
<div class="codehilite"><pre><span></span><code><span class="n">Value</span> <span class="nf">compute</span><span class="p">(</span><span class="k">const</span> <span class="n">Event</span> <span class="o">&amp;</span> <span class="n">ev0</span><span class="p">,</span> <span class="k">const</span> <span class="n">Event</span> <span class="o">&amp;</span> <span class="n">ev1</span><span class="p">);</span>
</code></pre></div>

<p>This version computes the distance between the two events as they are, without any preprocessing.</p>
<p><strong>Arguments</strong></p>
<ul>
<li><strong>ev0</strong><ul>
<li>The first event.</li>
</ul>
</li>
<li><strong>ev1</strong><ul>
<li>The second event.</li>
</ul>
</li>
</ul>
<p><strong>Returns</strong></p>
<ul>
<li>The Wasserstein distance between the two events.</li>
</ul>
<h4 id="ground_dists">ground_dists</h4>
<div class="codehilite"><pre><span></span><code><span class="n">std</span><span class="o">::</span><span class="n">vector</span><span class="o">&lt;</span><span class="n">Value</span><span class="o">&gt;</span> <span class="o">&amp;</span> <span class="n">ground_dists</span><span class="p">();</span>
<span class="k">const</span> <span class="n">std</span><span class="o">::</span><span class="n">vector</span><span class="o">&lt;</span><span class="n">Value</span><span class="o">&gt;</span> <span class="o">&amp;</span> <span class="n">ground_dists</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
</code></pre></div>

<p>To provide external ground distances, one should use this method to access a reference to the vector of distances used by the <code>NetworkSimplex</code> object. The user is free to resize this vector and copy/assign distances into it. Note that if <code>norm=false</code>, then  one of the events will generically get an extra particle to account for the difference in weights, and the external dists should account for this. The distances are expected to be provided in row-major (C-style) order.</p>
<p><strong>Returns</strong>
    - A (const) reference to the vector of ground distances used by the <code>NetworkSimplex</code> object.</p>
<hr />
<h3 id="access-results_1">Access Results</h3>
<p>Many of the methods share the same names as the <a href="#access-results">Python ones</a>. These are:</p>
<div class="codehilite"><pre><span></span><code><span class="n">Value</span> <span class="nf">emd</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">wasserstein</span><span class="o">::</span><span class="n">EMDStatus</span> <span class="nf">status</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>

<span class="n">std</span><span class="o">::</span><span class="n">vector</span><span class="o">&lt;</span><span class="n">Value</span><span class="o">&gt;</span> <span class="n">dists</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">std</span><span class="o">::</span><span class="n">vector</span><span class="o">&lt;</span><span class="n">Value</span><span class="o">&gt;</span> <span class="n">flows</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">Value</span> <span class="nf">flow</span><span class="p">(</span><span class="n">std</span><span class="o">::</span><span class="kt">ptrdiff_t</span> <span class="n">i</span><span class="p">,</span> <span class="n">std</span><span class="o">::</span><span class="kt">ptrdiff_t</span> <span class="n">j</span><span class="p">)</span> <span class="k">const</span><span class="p">;</span>
<span class="n">Value</span> <span class="nf">flow</span><span class="p">(</span><span class="n">std</span><span class="o">::</span><span class="kt">size_t</span> <span class="n">ind</span><span class="p">)</span> <span class="k">const</span><span class="p">;</span>

<span class="n">std</span><span class="o">::</span><span class="kt">ptrdiff_t</span> <span class="nf">n0</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">std</span><span class="o">::</span><span class="kt">ptrdiff_t</span> <span class="nf">n1</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">wasserstein</span><span class="o">::</span><span class="n">ExtraParticle</span> <span class="nf">extra</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">Value</span> <span class="nf">weightdiff</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">Value</span> <span class="nf">scale</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>

<span class="kt">double</span> <span class="nf">duration</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">std</span><span class="o">::</span><span class="kt">size_t</span> <span class="nf">n_iter</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>

<span class="n">std</span><span class="o">::</span><span class="n">pair</span><span class="o">&lt;</span><span class="n">std</span><span class="o">::</span><span class="n">vector</span><span class="o">&lt;</span><span class="n">Value</span><span class="o">&gt;</span><span class="p">,</span> <span class="n">std</span><span class="o">::</span><span class="n">vector</span><span class="o">&lt;</span><span class="n">Value</span><span class="o">&gt;&gt;</span> <span class="n">node_potentials</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
</code></pre></div>

<p><code>EMDStatus</code> and <code>ExtraParticle</code> are enums; see the <a href="../utils/#c">C++ Utils</a></p>
<hr />
<h3 id="getset-options_1">Get/Set Options</h3>
<p>These methods can be used to get the current settings of the <code>EMD</code> object or to set new ones. The getter methods are:</p>
<div class="codehilite"><pre><span></span><code><span class="n">Value</span> <span class="nf">R</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="n">Value</span> <span class="nf">beta</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="kt">bool</span> <span class="nf">norm</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="kt">bool</span> <span class="nf">do_timing</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="kt">bool</span> <span class="nf">external_dists</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="k">const</span> <span class="n">NetworkSimplex</span> <span class="o">&amp;</span> <span class="nf">network_simplex</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
<span class="k">const</span> <span class="n">PairwiseDistance</span> <span class="o">&amp;</span> <span class="nf">pairwise_distance</span><span class="p">()</span> <span class="k">const</span><span class="p">;</span>
</code></pre></div>

<p><code>external_dists()</code> indicates whether or not to use the internal <code>PairwiseDistance</code> object or to assume that the ground distances have been externally set. <code>network_simplex()</code> and <code>pairwise_distance()</code> access the underlying objects used to do the heavy lifting of the computation.</p>
<p>The setter methods are:</p>
<div class="codehilite"><pre><span></span><code><span class="kt">void</span> <span class="nf">set_R</span><span class="p">(</span><span class="n">Value</span> <span class="n">R</span><span class="p">);</span>
<span class="kt">void</span> <span class="nf">set_beta</span><span class="p">(</span><span class="n">Value</span> <span class="n">beta</span><span class="p">);</span>
<span class="kt">void</span> <span class="nf">set_norm</span><span class="p">(</span><span class="kt">bool</span> <span class="n">norm</span><span class="p">);</span>
<span class="kt">void</span> <span class="nf">set_do_timing</span><span class="p">(</span><span class="kt">bool</span> <span class="n">do_timing</span><span class="p">);</span>
<span class="kt">void</span> <span class="nf">set_external_dists</span><span class="p">(</span><span class="kt">bool</span> <span class="n">exdists</span><span class="p">);</span>
<span class="kt">void</span> <span class="nf">set_network_simplex_params</span><span class="p">(</span><span class="kt">unsigned</span> <span class="n">n_iter_max</span> <span class="o">=</span> <span class="mi">100000</span><span class="p">,</span>
                                <span class="n">Value</span> <span class="n">epsilon_large_factor</span> <span class="o">=</span> <span class="mi">1000</span><span class="p">,</span>
                                <span class="n">Value</span> <span class="n">epsilon_small_factor</span> <span class="o">=</span> <span class="mi">1</span><span class="p">);</span>
</code></pre></div>

<hr />
<h3 id="other-methods_1">Other Methods</h3>
<h4 id="preprocess">preprocess</h4>
<div class="codehilite"><pre><span></span><code><span class="k">template</span><span class="o">&lt;</span><span class="k">template</span><span class="o">&lt;</span><span class="k">class</span><span class="o">&gt;</span> <span class="k">class</span> <span class="nc">P</span><span class="p">,</span> <span class="k">typename</span><span class="p">...</span> <span class="n">Args</span><span class="o">&gt;</span>
<span class="n">EMD</span> <span class="o">&amp;</span> <span class="n">preprocess</span><span class="p">(</span><span class="n">Args</span> <span class="o">&amp;&amp;</span> <span class="p">...</span> <span class="n">args</span><span class="p">)</span>
</code></pre></div>

<p>Adds a preprocessor to the internal list. Each event will be preprocessed by the preprocessors in the order they were given. Currently, there is one preprocessor as part of Wasserstein, <a href="../utils/#centerweightedcentroid"><code>CenterWeightedCentroid</code></a>. Since this preprocessor takes no arguments, it can be added as:</p>
<div class="codehilite"><pre><span></span><code><span class="n">emd_obj</span><span class="p">.</span><span class="n">preprocess</span><span class="o">&lt;</span><span class="n">wasserstein</span><span class="o">::</span><span class="n">CenterWeightedCentroid</span><span class="o">&gt;</span><span class="p">()</span>
</code></pre></div>

<p>In general, any arguments to the preprocessor class are given as arguments to this method. A reference to the <code>EMD</code> object is returned.</p>
<h4 id="description_1">description</h4>
<div class="codehilite"><pre><span></span><code><span class="n">std</span><span class="o">::</span><span class="n">string</span> <span class="n">description</span><span class="p">(</span><span class="kt">bool</span> <span class="n">write_preprocessors</span> <span class="o">=</span> <span class="nb">true</span><span class="p">)</span>
</code></pre></div>

<p>Returns a string that describes the <code>EMD</code> object.</p>
<p><strong>Arguments</strong></p>
<ul>
<li><strong>write_preprocessors</strong><ul>
<li>Whether or not to include preprocessors in the description. There are currently no preprocessors included in the Wasserstein package but this may change in the future.</li>
</ul>
</li>
</ul>
<p><strong>Returns</strong></p>
<ul>
<li>The description of the <code>EMD</code> object.</li>
</ul>
<h4 id="clear_1">clear</h4>
<div class="codehilite"><pre><span></span><code><span class="kt">void</span> <span class="n">clear</span><span class="p">()</span>
</code></pre></div>

<p>Frees some memory in use by the <code>EMD</code> object. This should not normally need to be called by the user.</p>

            </div>
          </div><footer>
    <div class="rst-footer-buttons" role="navigation" aria-label="Footer Navigation">
        <a href="../../releases/" class="btn btn-neutral float-left" title="Release Notes"><span class="icon icon-circle-arrow-left"></span> Previous</a>
        <a href="../emds/" class="btn btn-neutral float-right" title="PairwiseEMD">Next <span class="icon icon-circle-arrow-right"></span></a>
    </div>

  <hr/>

  <div role="contentinfo">
    <!-- Copyright etc -->
      <p>Copyright (C) 2019-2021 Patrick T. Komiske III</p>
  </div>

  Built with <a href="https://www.mkdocs.org/">MkDocs</a> using a <a href="https://github.com/readthedocs/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>.
</footer>

        </div>
      </div>

    </section>

  </div>

  <div class="rst-versions" role="note" aria-label="Versions">
  <span class="rst-current-version" data-toggle="rst-current-version">

        <span>
          <a href="https://github.com/thaler-lab/Wasserstein/" class="fa fa-github" style="color: #fcfcfc"> GitHub</a>
        </span>


      <span><a href="../../releases/" style="color: #fcfcfc">&laquo; Previous</a></span>


      <span><a href="../emds/" style="color: #fcfcfc">Next &raquo;</a></span>

  </span>
</div>
    <script>var base_url = '../..';</script>
    <script src="../../js/theme_extra.js" defer></script>
    <script src="../../js/theme.js" defer></script>
      <script src="../../js/extra.js" defer></script>
      <script src="https://polyfill.io/v3/polyfill.min.js?features=es6" defer></script>
      <script src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-chtml.js" defer></script>
      <script src="../../search/main.js" defer></script>
    <script defer>
        window.onload = function () {
            SphinxRtdTheme.Navigation.enable(true);
        };
    </script>

</body>
</html>
